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We determine the energetically lowest lying states in the BEC-BCS crossover regime of s-wave 
interacting two-component Fermi gases under harmonic confinement by solving the many-body 
Schrodinger equation using two distinct approaches. Essentially exact basis set expansion techniques 
are applied to determine the energy spectrum of systems with A*' = 4 fermions. Fixed-node diffusion 
Monte Carlo methods are applied to systems with up to A'' = 20 fermions, and a discussion of 
different guiding functions used in the Monte Carlo approach to impose the proper symmetry of the 
fermionic system is presented. The energies are calculated as a function of the s-wave scattering 
length as for A*' = 2 — 20 fermions and different mass ratios k of the two species. On the BEC and 
BCS sides, our energies agree with analytically-determined first-order correction terms. We extract 
the scattering length and the effective range of the dimer-dimer system up to k = 20. Our energies 
for the strongly-interacting trapped system in the unitarity regime show no shell structure, and are 
well described by a simple expression, whose functional form can be derived using the local density 
approximation, with one or two parameters. The universal parameter ^ for the trapped system for 
various k is determined, and comparisons with results for the homogeneous system are presented. 



PACS numbers: 

I. INTRODUCTION 

Advances in trapping and cooling have spawned the 
experimental realization of ultracold externally-confined 
two-component Fermi gases with controllable interac- 
tion strengths. Using these impurity-free systems, the 
crossover from a weakly-attractive atomic Fermi gas 
through a strongly-interacting unitarity regime to a 
weakly-repulsive molecular Bose gas has been investi- 
gated i^, [3), 1J] . Our increased understanding of these 
systems relates to the study of neutron matter and po- 
tentially that of high-Tc superconductors, in addition to 
the field of ultracold atomic gases. All of these systems 
are controlled by similar pairing mechanisms, although 
at much different densities. 

To date, experimental studies of the BEC-BCS 
crossover with ultracold atomic gases have been re- 
stricted to fermions in different hyperfine substates. In 
this case, the "spin-up" and "spin-down" fermions have 
equal masses and experience equal trapping frequencies. 
Currently, the simultaneous trapping and cooling of dif- 
ferent atomic fermionic species is being pursued in a num- 
ber of laboratories. This motivates us to investigate how 
the BEC-BCS crossover physics changes with the mass 
ratio K of the two atomic species. Our goal is to develop a 
microscopic understanding of these intricate many-body 
systems. To this end, we consider trapped systems with 
varying number of particles N , and relate them to the 
homogeneous systems through the local density approx- 
imation (LDA). This illuminates the transition from the 
few-body to the many-body physics of an ultracold Fermi 
gas. 

For a given short-range two-body potential, the 
stationary Schrodinger equation for a trapped two- 
component Fermi gas has a rich eigenspectrum, which in 



some cases includes deeply- and/or weakly-bound clus- 
ter states as well as "ground" and highly-excited gas- 
like states. In general, the eigenstates of two-component 
Fermi gases with short-range interactions can be sepa- 
rated into two classes: universal states that do not (or 
only weakly) depend on the details of the two-body po- 
tential [1] , and non-universal states that depend notably 
on the details of the two-body potential. The eigenstates 
of the four-fcrmion system with equal masses, e.g., fall 
into the former class, provided the range of the two-body 
potential is sufficiently small; in this case, the proper- 
ties of the system are to a very good approximation de- 
termined by a single parameter, the s-wave scattering 
length as- For large mass ratios, however, non- universal 
bound trimer states exist @, 0, [Ml . A description of these 
states requires a three-body parameter, which depends 
on the short-range physics. In some cases, non-universal 
bound clusters consisting of four or more fermions may 
exist. In this work, we do not analyze the properties of 
such non-universal states but instead study the proper- 
ties of states that depend at most weakly on the short- 
range physics. In particular, we determine the BEC-BCS 
energy crossover curve, which is defined in Sec. IIII Al 
by solving the stationary Schrodinger equation for var- 
ious mass ratios n. An analysis of the stability of two- 
component Fermi systems with large mass ratios, includ- 
ing molecular Bose gases created from two-component 
Fermi gases, is beyond the scope of this paper. 

Consider the stationary solutions of the four-particle 
system as a function of n in the BEC-BCS crossover. 
The Schrodinger equation is solved in two distinct ap- 
proaches: a basis set expansion technique that utilizes 
correlated Gaussians (CG) and a fixed-node diffusion 
Monte Carlo (FN-DMC) approach. The dimer-dimer 
scattering length add and the dimer-dimer effective range 
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rdd emerge as a function of k. A surprisingly large rdd 
is found, which is likely to be an important input pa- 
rameter in the BEC many-body theory. Furthermore, a 
detailed comparison of the results throughout the entire 
crossover regime permits a non-trivial test of the nodal 
surface employed in the FN-DMC approach, and it con- 
veys information about the symmetry of the many-body 
wave function. Extension of our FN-DMC calculations to 
larger numbers of particles also probes the validity range 
of the analytically-determined limiting behaviors in the 
deep BCS and BEC regimes. In the strongly-interacting 
unitarity regime, the LDA relates the trapped system 
properties to those of the homogeneous system. Finally, 
our FN-DMC energies should allow for a stringent test 
of numerically less involved approaches such as density 
functional treatments [91]. 

Section [IT] introduces the Hamiltonian of the trapped 
Fermi system, and the numerical approaches ap- 
plied to solve the corresponding stationary many-body 
Schrodinger equation. Section IIIII presents our results 
for the energetics and the interpretation of the results of 
weakly- and strongly-interacting Fermi systems with up 
to 20 atoms. Finally, Sec. IIVI concludes. 



the BEC-BCS crossover using magnetic field strengths 
for which the p-wave interactions are non-resonant. 

The Hamiltonian given by Eq. ([1]) is characterized by 
four different length scales: the range i?o of the interac- 
tion potential, the s-wave scattering length a^, and the 
two oscillator lengths a'H — ^Jh/{mjU)j), j = I and 2. 
Throughout, we are interested in the regime where Rq is 
much smaller than the oscillator lengths a|,"2i or equiva- 
lently, where the system is dilute with respect to Rq, i.e., 
n(0)i?Q ^ 1, where n(0) denotes the peak density. In 
this regime, the properties of the universal states are ex- 
pected to be independent of the details of the two-body 
potential. For a given as, we numerically test whether a 
state is universal by calculating its energy for various Rq- 
The condition Rq ^ a^f^^ implies that the numerical ap- 
proaches chosen for solving the many-body Schrodinger 
equation have to be able to govern the physics occuring 
at at least two different length scales. As we illustrate 
below, the CG and FN-DMC approaches are able to do 
so. 



B. Correlated Gaussian approach 



II. HAMILTONIAN AND NUMERICAL 
APPROACH 

A. Hamiltonian 

For N harmonically-trapped Fermi atoms divided 
equally into two species, the Hamiltonian is given by 
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Here, unprimed indices label mass mi and primed indices 
mass 7712 fermions, and N is assumed to be even. The 
mass ratio k is defined by mi /m2 and throughout we take 
TOi > m2. In Eq. ((T|), uji and UJ2 denote angular trapping 
frequencies, and ri the position vector of the ith fermion. 

We adopt two purely attractive short-range model 
potentials for the interaction between unlike fermions: 
a Gaussian interaction potential V{r), V{r) — 
— Vb exp(— r^/(2i?§)), and a square well interaction po- 
tential V{r), V{r) — ~Vo for r < Rq and otherwise. 
For a fixed range Rq of the two-body potential V{r), the 
depth Vq, Vq > 0, is adjusted until the s-wave scattering 
length Os assumes the desired value. For negative (or pos- 
itive) as , Vo and Rq are chosen so the potential supports 
no (or one) two-body s-wave bound state. Throughout 
this paper, we treat the like atoms as non-interacting. 
This is justified because the interactions between like 
atoms are only non- negligible very close to a p-wave Fes- 
hback resonance. All experiments to date have studied 



The correlated Gaussian methodflOl. [Tl| is a powerful 
tool to study few-body systems. Recently, the CG ap- 
proach has been applied to the four-fermion system with 
equal masses |12i |. Here, we analyze the properties of this 
system from a somewhat different point of view and ad- 
ditionally consider the unequal mass system. As in the 
previous work, we treat equal trapping frequencies, i.e., 
LUi — 072, so that the center-of-mass motion separates. 
To further reduce the dimensionality of the problem, we 
restrict ourselves to states with vanishing total angular 
momentum L and positive parity P. We expand the 
= 0"^ states in terms of correlated Gaussian basis 
functions ^j-, which depend on the six interparticle dis- 
tances and the center of mass vector. 
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Here, the Cj- denote expansion coefficients. Each basis 
function <I>j- is written as a product of the ground state 
center-of-mass wavefunction and a symmetrized product 
of Gaussian functions for each interparticle distance vec- 
tor [12] ■ The widths of these Gaussians can be different 
for each interparticle distance, giving us six parameters 
for each basis function. These parameters are in Eq. ^ 
collectively denoted by d. The simple functional form of 
the wave function ip, Eq. allows the analytical deter- 
mination of all matrix elements if the two-body interac- 
tion potential is taken to be a Gaussian (see Sec. Ill Ap . 

The parameter vector d that characterizes each basis 
function are selected semi-randomly. Typically, the 

components of d vary from a fraction of the range Rq to 
a few times the interparticle distance in the noninteract- 
ing limit. The basis functions can be roughly separated 
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into three types. The first type has all the components 
of d of the order of the trap lengths and is suitable to 
describe gas-like states. The second type has one or two 
small d components while the others take values of the 
order of the trap lengths; these basis functions carry a 
large weight when describes states that consist — to a 
good approximation — of two bound dimers or a dimer 
and two free atoms. The third type has more than two 
small d components, and is suitable to describe compara- 
tively tightly-bound three- and four-body states. In gen- 
eral, all three types of basis functions are needed to ac- 
curately describe the entire = 0+ spectrum of the 
four-fermion system. For equal masses, however, we find 
that the third type carries negligible weight, owing to the 
absence of molecular three- or four-body states. As the 
mass ratio increases, more basis functions of the third 
type need to be included. We carefully check the con- 
vergence of the energies by varying the total number of 
basis functions used, and by varying the basis functions 
included in the expansion. We find that of the order 
of 10'* basis functions suffice to accurately describe the 
eigenfunctions of interest. 

The basis functions introduced above are not linearly 
independent. To eliminate the linear dependence in our 
basis set, we diagonalize the overlap matrix and elimi- 
nate the basis functions with the lowest eigenvalues up 
to a certain cutoff. The remaining basis functions are 
then used to construct a new orthogonal basis set. Fi- 
nally, the eigenspectrum is obtained by diagonalizing the 
corresponding Hamiltonian matrix. 



C. Fixed- node diffusion Monte Carlo approach 

To treat up to = 20 fermions, solutions of the 
Schrodinger equation are determined by the FN-DMC 
method [13l . Il4| . In this method, the proper fermionic 
antisymmetry is imposed through the use of a so-called 
guiding function ipx, which depends on the coordinates of 
all particles. To within statistical uncertainties, the FN- 
DMC algorithm provides an upper bound to the exact 
ground state energy, i.e., to the lowest-lying state with 
the same symmetry as ipT- Note that all our FN-DMC 
calculations are performed for the square well interaction 
potential. While there is no technical problem in extend- 
ing the FN-DMC calculations to the Gaussian interaction 
potential, the guiding functions ipT are most readily de- 
termined and evaluated for the square well potential. 

Two different guiding functions ipr considered in this 
work are: 

N/2 N/2 
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TABLE I: FN-DMC energies E in units of hio as a function 
of N, N — 2 — 20, for the two-component Fermi system with 
a^^^ = a^f^^ at unitarity for k, — 1 and 8. Statistical uncertain- 
ties of the energies are reported in round brackets. Note that 
no nodal approximation needs to be made for N = 2. 
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and 

N/2 

^pT2 = *JV/ (ri , • • • , rN/2' ) X Y[ f{ni' ) , (4) 

Here $ denotes the ground state harmonic oscillator or- 
bital, A is the antisymmetrizer, and 5* at/ denotes the 
wave function of A^ trapped non-interacting fermions. 
Following Ref. 15], the pair function / is constructed 
from the free-space zero-energy scattering and the free- 
space bound state solution of the two-body square-well 
interaction potential for negative and positive s-wave 
scattering length a^, respectively. In Eq. / coincides 
with / for small r and is matched smoothly to a non- 
zero constant at larger r. This matching to a non-zero 
constant ensures that the product over all pair functions 
/ is always non-zero. Thus, the nodal structure of iIjt2 
coincides with that of the non-interacting Fermi gas. In 
contrast, the nodal surface of ipTi is constructed by an- 
tisymmetrizing a product of pair functions p^ . 

To assess the accuracy of our MC code, we determine 
the energy of the two-body system with lui = uj2 and 
TOi = 1712 semi-analytically. We separate off the center- 
of-mass motion, and write the eigenfunctions of the 
Schrodinger equation for the relative coordinate in terms 
of hypergeometric functions. The resulting eigenequa- 
tion results in an energy of E — 2.00200?ia). Since the 
two-body wave function is nodeless, the DMC energy for 
N = 2 (see Table HI is expected to be exact. Indeed, 
the DMC energy agrees to within the statistical uncer- 
tainty with the energy determined semi-analytically. A 
detailed comparison of the FN-DMC and CG energies for 
the four-fermion system, which allows the quality of the 
nodal surface employed in the FN-DMC calculations to 
be assessed, is presented in Sec. IIIIBl 

In the non-interacting case, i.e., for as = 0, the guid- 
ing function '0^2 with f(r) = 1 coincides with the exact 
eigen function. For weakly-attractive Fermi systems, the 
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attractive nature of the two-body potential introduces 
correlations but does, to a good approximation, leave 
the nodal surface unchanged. Indeed, we find that the 
variational energy for in this regime is nearly in- 
distinguishable from the FN-DMC energies, indicating 
that the Jastrow product over all pair functions accounts 
properly for the two-body correlations of the system and 
that three- and higher-order correlations are negligible. 

For small positive Os, on the other hand, compara- 
tively strongly-bound two-body dimers exist and the sys- 
tem is expected to form a molecular Bose gas of dimers. 
Such a system is not even qualitatively described cor- 
rectly by the guiding function ipTi, which assumes that 
every spin-up fermion is "simultaneously" correlated with 
every spin-down fermion [l6l | . The guiding function , 
instead, is much better suited to describe a Fermi gas 
that behaves as a weakly-interacting molecular Bose gas. 
iJjt2 correlates the first spin-up fermion with the first 
spin-down fermion, the second spin-up fermion with the 
second spin-down fermion, and so on, and then anti- 
symmetrizes this "paired state". The guiding function 
ipT2 is expected to accurately describe the system when 
the size of the dimer pairs becomes small compared to 
the oscillator lengths. 

Finally, in the strongly- interacting regime, i.e., for 
jflsl — > oo, it is not a priori clear which of the two guid- 
ing functions provides a better description of the system. 
Section IIIII discusses this in more detail, and also com- 
ments on additional aspects of the choice of the guiding 
functions related to the existence of non-universal trimer 
states. 



III. RESULTS 

A. Energy crossover curve: Definition and general 
consider at ions 

Throughout this work, we are interested in describ- 
ing the crossover from a weakly-repulsive to a weakly- 
attractive trapped two-component Fermi system with N 
particles. This BEC-BCS crossover can be characterized 
by the normalized energy crossover curve A^'' , 

_ EiN)-NE{2)/2 

which depends on a^, k and N. In Eq. ([5]), E{N) denotes 
the energy of the A^-fermion system, ui = {uji + a-'2)/2 is 
the average frequency and A is defined through the energy 
Eni of N non-interacting fermions, 

Eni = iX + SN/2)hLd. (6) 

The values of A for the first few closed-shell systems are 
listed in the second column of Table [Til ^iv^ equals one 
on the deep BCS side (small \ag\ and < 0), and zero 
on the deep BEC side (small Ug and > 0). Since the 
energy E{2) of N/2 trapped dimer pairs is subtracted 



from the total energy E{N) of the system, the energy 
crossover curve A^"*, Eq. ([5]), is expected to be indepen- 
dent of the details of the two-body potential if the range 
Ro is much smaller than the average interparticle spacing. 
The energy crossover curve defined here for the trapped 
system is the analog of the BEC-BCS crossover curve of 
the homogeneous system (see, e.g., Refs. [13, [3 [3 for 
pioneering work based on the mean-field BCS equations, 
and Figs. 1 and 2 of Ref. [l^ for a determination of the 
crossover curve for the homogeneous system by the FN- 
DMC method). 

As indicated above, A^^ depends on the scattering 
length as, the number of particles N, the ratio between 
the two masses, and the ratio between the two frequen- 
cies. Thus, an exhaustive study of the whole parameter 
space of trapped two-component Fermi systems by first 
principle methods is impossible. This paper considers 
two different scenarios: i) the trapping frequencies are 
set to coincide, i.e., lui = lo2, while k, N and as are 
varied (see Sees. IIII Bl and IIII C|) . and ii) the oscillator 

lengths are set to coincide, i.e., a'f^^ = a^j^^ , and l/\as\ is 
set to 0, while N and k are varied (see Sec. IIII D|) . 

Our motivation for considering scenario i) is as fol- 
lows. In the deep BEC regime, the fermionic system is 
expected to form a molecular Bose gas whose behaviors 
are to a good approximation determined by the dimer- 
dimer scattering length add- The dimer-dimer scatter- 
ing length can be extracted quite readily for different 
K from our four-body energies, provided the center-of- 
mass motion decouples (see Sec. IIIIBj) . For unequal fre- 
quencies, the center-of-mass motion does not decouple 
and the extraction of the dimer-dimer scattering length 
would be more involved. Section [III CI extends the study 
of the four-fermion system with equal frequencies to sys- 
tems with more particles to illustrate that the behaviors 
of the larger systems in the deep BEC regime are also 
governed to a good approximation by the dimer-dimer 
scattering length. While this has been shown to be the 
case previously for the homogeneous system with equal 
masses [15| , our calculations illustrate that — as might be 
expected — the many-body physics of unequal mass sys- 
tems in the deep BEC regime is also to a good approxima- 
tion governed by a single few-body parameter, the dimer- 
dimer scattering length. In the more strongly-interacting 
regime, our equal frequency calculations for k > 1 pro- 
vide insights into the behaviors of systems whose den- 
sities are not fully overlapping. Thus, mass-imbalanced 
Fermi system may behave in certain respects similar to 
population-imbalanced Fermi systems. 

Our primary motivation for considering scenario ii) is 
to connect the behaviors of the trapped system with those 
of the homogeneous system using the LDA. The energy of 
the homogeneous system at unitarity is related to the en- 
ergy of the non-interacting system by a universal param- 
eter ^. By calculating the energies of the trapped system 
at unitarity for equal frequencies and equal masses, we 
quantify how well the LDA describes small trapped sys- 
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FIG. 1: (Color online) Energy crossover curve A*"' for 

uji — L02 as a function of a!'^^'' /as for (a) k = 1 and (b) 
K = 8. Solid lines are calculated by the CG approach, and 
circles and crosses by the FN-DMC method using ijjTi and 
V'T2, respectively. 



terns. Unpublished results for the homogeneous system 
with equal spin-up and spin-down densities but unequal 
masses indicate that the universal parameter ^ depends 
weakly on the mass ratio k (20| . These unequal mass 
results for the homogeneous system cannot be straight- 
forwardly connected to those of the trapped system since 
the densities of the two trapped species do not necessarily 
fully overlap. In the non-interacting limit, the densities 
of the two species coincide when the trapping lengths are 
equal, i.e., a'^l = o'^l- At unitarity, we find that unequal 

mass systems with a'j^l — aj^^ have a small, though non- 
negligible, density mismatch. Despite this small density 
mismatch, we apply the LDA a"j^^ = aJ^^J and relate the 
universal parameter of the homogeneous system to the 
energetics of the trapped system. 



B. Energy crossover curve for A'^ = 4 

This subsection presents our results for the energy 
crossover curve for four trapped fermions calculated by 
the CG and FN-DMC approaches for equal frequencies, 
i.e., uji ~ LJ2 = iij, and varying mass ratio k. 

Figure [1] shows the energy crossover curve for 

four fermions as a function of cLhlf' /ds calculated by the 
CG and FN-DMC approaches for (a) k = 1 and (b) 
K = 8. Here, the oscillator length a^j^l^^ is defined in 
terms of the reduced mass /i — 17111712/ (rni + TO2), i.e., 

^ho^ ~ VWi^i^- The solid lines in Fig. [T] are ob- 
tained using E{4) calculated by the CG approach, while 
circles and crosses are obtained using E{i) calculated 
by the FN-DMC approach using ipTi and ipT2, respec- 
tively. The ranges Rq of the two-body potentials used in 
Fig. [T]are much smaller than the oscillator lengths, i.e.. 



Rq ~ O.Ola\'^^\ From our CG energies for different Rq, 

we estimate that the Aj^' shown in Fig. [1] deviate by at 
most 1% from the corresponding curves for zero-range 
interactions. For mi — 7712, e.g., the energy at unitarity 
calculated by the CG approach for the Gaussian inter- 



and 



action potential is E — 5.027fkj for Rq = O.Ola^^ 

E = 5.099to for Rq ~ COSalj^^^-* [5]. For comparison, 
the FN-DMC energy for the square well potential with 
i?o = O.Ol^o''^ is £; = 5.069(9) (see Table H]), which is in 
good agreement with the energies calculated by the CG 
approach. 

As expected, the energy crossover curve connects the 
limiting values of one on the BCS side and zero on the 
BEG side smoothly. Importantly, the lowest FN-DMC 
energies and the CG energies agree well, which implies 
that the functional forms of t/^^i and V'T2 are adequate. 
For equal masses [panel (a)], the FN-DMC energies at 
unitarity calculated using the two different agree ap- 
proximately. For unequal masses [panel (b)], in contrast, 
the nodal surface of "072 leads to a lower energy at uni- 
tarity than that of '0ri , and the crossing point between 
the energies calculated using and 1/^2 moves to the 
BEG side. This can be understood by realizing that the 
densities of the heavy and light particles do not overlap 
fully, leading to a reduced pairing. 

The CG approach in our current implementation (see 
Sec. Ill Bl) allows for the determination of the complete 
= 0"*" energy spectrum. If we use short-range Gaus- 
sian two-body potentials that support no two-body s- 
wave bound state for negative Og and one two-body s- 
wave bound state for positive a^, the four-body energy 
that enters the calculation of the energy crossover curves 
shown in Fig. [T] is the true ground state of the sys- 
tem, i.e., no energetically lower-lying bound trimer or 
tetramer states with = 0"*" symmetry exist. For 
larger mass ratios, bound trimer states exist. The mass 
ratio at which these non-universal trimer states appear 
depends on the range i?o of the two-body potential em- 
ployed. In the regime where three-body bound states 
exist, the four-body spectrum calculated by the CG ap- 
proach contains also universal states which are separated 
by approximately 2to and which can be best described 
as two weakly-interacting composite bosons. For fixed 
a Si a,s > 0, the energy of these "dimer-dimer states" 
changes smoothly as a function of k even in the regime 
where bound trimer states appear. In the following, we 
use these dimer-dimer states to extract the dimer-dimer 
scattering length as a function of k up to k = 20. 

When as is small, the four fermions form two bosonic 
molecules of mass M, where M — nii + TO2, which inter- 
act through an effective molecule-molecule potential [8J. 
To model this effective dimer-dimer potential (the ex- 
act functional form is unknown) , we introduce a regular- 
ized zero-range potential V{r) [22|, V{r) = gS{r){d/dr)r, 
whose scattering strength g is parameterized by the scat- 
tering length add and the effective range rdd of the dimer- 
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FIG. 2: (Color online) Four-body energies of the three en- 
ergetically lowest-lying dimer-dimer states as a function of 



fls/aLt^' for K 



Panel (a) shows the energetically lowest 



lying energy level (i = 0), panel (b) the energetically second 
lowest {i = 1) and panel (c) the energetically third lowest 
state (i = 2). Circles and crosses show our CG and FN-DMC 
results, respectively. Solid lines show the zero-range model 
results. 



dimer system, i.e., 

4:TTh^ add 
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M 



MEtbTddadd 
2f9 



(7) 



Here, Etb denotes the relative energy of the two-boson 
system, i.e., the total energy with the center-of-mass 
contribution subtracted. It has been shown previ- 
ously [13, nil that the inclusion of the energy depen- 
dence of the scattering length notably extends the va- 
lidity regime of the zero-range pseudopotential when ap- 
plied to describe the scattering of two atoms under exter- 
nal confinement. By comparing the "dimer-dimer energy 
levels" of the four-fermion system with the energies of 
two mass M bosons in a trap interacting through this 
energy-dependent zero-range potential [22, 23, l2J|, we 
determine add and Vdd- 

To illustrate this procedure, circles in Figs. O^a) 
through (c) show the three energetically lowest-lying 
dimer-dimer energy levels, referred to as £'j(4) with i = Q 
through 2, with the center-of-mass energy and the dimer- 
binding energy subtracted for k = 8 obtained by the CG 
approach as a function of . Solid lines show the energy 
levels obtained by fitting these four-body energies by the 
two-boson energies obtained using the energy-dependent 
zero-range pseudopotential {add and rdd are treated as 
fitting parameters). We find that inclusion of the effec- 
tive range rdd extends the validity regime over which the 
four-fermion system can be described by the two-boson 
model and additionally allows for a more reliable deter- 
mination of a (id- Figure [2] illustrates that the two-boson 
spectrum reproduces the dimer-dimer states of the four- 




FIG. 3: (Color online) Circles and crosses show add/ as as a 
function of k, extracted from the four-fermion CG and FN- 
DMC energies, respectively. For comparison, a solid line 
shows the results from Fig. 3 of Ref. @]. Diamonds and 
squares show Vdd/as extracted from the four-fermion CG and 
FN-DMC energies, respectively. 



fermion spectrum well over a fairly large range of atom- 
atom scattering lengths a^. For comparison, crosses in 
Fig. [2ja) show the corresponding FN-DMC energies for 
the energetically lowest-lying dimer-dimer state. We did 
not attempt to construct a guiding function that would 
allow for the determination of excited dimer-dimer states. 
We find that the FN-DMC energies are slightly larger 
than the CG energies and that the deviation increases 
with increasing as- Presumably, this can be attributed to 
the functional form of the nodal surface used in the FN- 
DMC calculations, which should be best in the very deep 
BEC regime. The increasing deviation between the FN- 
DMC and CG energies with increasing explains why 
the effective range predicted by the analysis of the FN- 
DMC energies is somewhat larger than that predicted by 
the analysis of the CG approach (see discussion of Fig. [3] 
below) . 

Circles and crosses in Fig. [3] show the resulting dimer- 
dimer scattering lengths add extracted from the energies 
calculated by the CG and the FN-DMC approach, re- 
spectively, as a function of n. For all mass ratios con- 
sidered in Fig. [3l we include up to three dimer-dimer 
energy levels in our analysis of the CG results, and only 
the lowest dimer-dimer level in our analysis of the FN- 
DMC results. Our dimer-dimer scattering lengths agree 
well with those calculated by Petrov et al. within a zero- 
range framework Q (solid line in Fig. [3]) . The calcula- 
tions by Petrov et al., performed for the free and not the 
trapped four-fermion system, terminate at k ~ 13.6, be- 
yond which a three-body parameter is needed to solve the 
four-body equations within the applied framework. Our 
calculations show the existence of deeply-bound "plung- 
ing" states, which consist of a trimer plus a free atom. 
This signals a qualitative change of the energy spectrum, 
in agreement with Petrov et al. 0. At the same time, 
our calculations for finite-range potentials predict that 
add continues to increase smoothly when the mass ratio 
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TABLE II: Values of A and for the four smallest closed- 
shell two-component Fermi systems with equal frequencies. 



iV 



A 



2 

8 

20 
40 





6 

30 
90 



2 (4+23k.+4k^) 

2 5(20+50k-I-249k^+50«:^-|-20«:'^) 
4(1-1- fc)* 

2 5(40+450K-l-306/t^+2795K^-f-306K^-l-450K^+40K°) 
4,(1+k)>> 



FIG. 4: (Color online) Energies for N = 8 and small \as\. (a) 

E{8)/Eni as a function of as/a^i^^^ for wi = tJ2 = 'i' for k = 1 
(circles and solid line) and n — 8 (crosses and dashed line). 
Symbols are calculated by the FN-DMC method using ^t2, 
and lines using the first order correction, Eq. (|10|l . (b) Ag"^' 
as a function of as/a)^P for uji — lj2 ~ for k = 1 (circles 
and solid line) and k = 8 (crosses and dashed line). Symbols 
are calculated by the FN-DMC method using 'ipri, and lines 
using the first order correction, Eq. 

K exceeds 13.6. This can possibly be explained by the fact 
that the presence of the external confining potential may 
"wash out" some of the features present in the free-space 
system. As already mentioned, the study of the stabil- 
ity of the four-fermion system, consisting of two dimcrs, 
with large mass ratios is beyond the scope of this work. 

Diamonds and squares in Fig. [3] show the effective 
range rdd extracted from our CG and FN-DMC ener- 
gies, respectively. We estimate the uncertainty of rdd 
obtained from the CG approach to be about 10%, and 
quite a bit larger for that extracted from the FN-DMC 
energies. Figure [3] shows that the ratio rdd/ add increases 
from about 0.2 for k = 1 to about 0.5 for k = 20. While 
earlier work already suggested that the dimer-dimer po- 
tential may be best characterized as a broad soft-core 
potential [8|, implying a non-negligible value for the ef- 
fective range rdd, our work makes the first quantitative 
predictions for rdd as a function of k. The large value 
of rdd suggests that effective range corrections may need 
to be considered in analyzing the physics of molecular 
Fermi gases. 



C. Weakly-interacting limits for > 4 

We now apply the FN-DMC method to larger systems 
with wi = ^2 = w, focussing on the deep BCS and BEG 
regimes where |as| is small. Figure HJa) shows the total 
energy E{N) for N = 8 particles, divided by the energy 
Ejfi of the non-interacting system, for k = 1 and 8 as a 
function of as/a-^o^ for small \as\. The FN-DMC ener- 
gies, calculated using ■!/'T2, are shown by symbols (circles 
for K = 1 and crosses for k = 8). 

For comparison, we calculate the energy of a weakly- 
attractive closed-shell Fermi system with equal frequen- 
cies ((D = wi = LJJ2) in first order perturbation theory. We 



assume that the unlike fermions are interacting through 
the Fermi pseudopotential Vs{r, r') = '^^^^^6{r-r') [2|. 
Applying perturbation theory to the non-interacting two- 
component Fermi gas with unequal masses but equal fre- 
quencies, the energy becomes E « E^i + E^^^, where the 
first order energy correction e'>^^ for closed-shell systems 
can be written as 

Eill = ^-^ [p„.Ar)PmAr)dr. (8) 

Here, Prm denotes the density of a single non-interacting 
mass rui component (i = 1 and 2). The integration in 
Eq. ([5]) can be performed analytically, resulting in a sim- 
ple expression for the first-order energy correction, 

Elll = f^C^N^, (9) 

where denotes a constant that depends on N and k. 
Altogether, we obtain 

ENi + hioC^^. (10) 

The values of for the first four closed-shell systems 
are summarized in the third column of Table [III For 
completeness, the second column summarizes the values 
of A that determine the energy E^i of the non-interacting 
system [see Eq. 

The first order correction, Eq. (jlOp . is shown by solid 
and dashed lines in Fig. [4lja) for k = 1 and 8, re- 
spectively; it describes the interacting system well for 
\as\ < 0.25a'^//\ or equivalently, for fc^a^ < 0.6. Here, 
kp denotes the wave vector at the trap center of a non- 
interacting system of N fermions with mass 2/i eval- 
uated within the Thomas Fermi approximation, kp = 

V2(37V)i/Vaio^^- Additional corrections can be derived 
within a renormalized scattering length framework (26j . 

We now turn to the small Os limit for N = 8. Circles 
(k = 1) and crosses (k — 8) in Fig.[51[b) show the energy 
crossover curve Ag'"'' in the BEG regime for Co = uji = uj2, 
where E{S) is calculated by the FN-DMC method, as 
a function of as/a)^P , where ajj^ — ^fij (Mw); in the 

small fls regime, a'^^ is the relevant characteristic oscil- 
lator length. Treating the 7V-fermion system as a bosonic 
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FIG. 5: (Color online) Circles and crosses show the FN-DMC 
energies E{N) in units of hCu as a function of at unitarity 
for K — 1 {uji =012) and k = 8 (w2 = /f^i), respectively. Solid 
lines show a fit of the FN-DMC energies to Eq. (fT3)l. 



gas consisting of N/2 mass M molecules and applying 
first order perturbation theory using a Fermi pseudopo- 
tential [2^ , the energy of the system with u) — uii — UJ2 
reads 



NiN - 2) /2 add 



AM) ■ 



(11) 



SoHd and dashed lines in Fig. [l{b) show Ag for k = 1 
and 8, respectively, calculated using Eq. (fTT|) . To plot the 
expansion, we use the dimer-dimer scattering length add 
calculated by the CG approach. For both mass ratios, the 
agreement between the FN-DMC energies and the first 
order correction is good for Os < O.ba^^P . Figure [31^ a) 
thus illustrates that the behavior of the Fermi system 
depends to a first approximation only on add if '^s/o.'hP is 
sufhciently small. Inclusion of effective range corrections 
may improve the agreement but is beyond the scope of 
this paper. 

We checked that the behaviors discussed here for N — 
8 also hold for = 20 particles. 



D. Energetics at unitarity 

This section considers the strongly-interacting unitary 
regime, where the atom-atom scattering length is infi- 
nite. To ensure large overlap of the densities of the two 
species, we choose the trapping frequencies uji and UJ2 so 



that a 



(1) 

ho 



,(2 



'J (see Sec. IIII K\ for a discussion). Circles 
and crosses in Fig. [5] show our FN-DMC energies E{N) 
at unitarity as a function of N for n — 1 and 8, respec- 
tively, while Table [J fists tfie FN-DMC energies. In these 
calculations, the range Rq of the square well potential 
used to describe the interaction between unlike fermions 
is set to i?o = O.Ola^^j. The energies for = 4 are cal- 
culated using tpT2', usage of tpTi leads to slightly higher 
energies. As discussed already in Sec. IIIIBl the four- 
body FN-DMC energy for equal masses agrees well with 



the corresponding CG energy. The energies for N > 6 
are calculated using the guiding function ipxi- For ex- 
ample, usage of the guiding function 'ipT2 gives an energy 
of f2.64(2)fia; for A^ = 8 and k — 1 (which is, taken 
the statistical errorbars into account, just slightly higher 
than the energy calculated using V'Ti; see Table |T]), and 
an energy of 43.2{l)hLu for A^ = 20 and k = 1 (which is 
notably higher than the energy calculated using ■f/'Ti; see 
Table III). 

For A^ > 8, our energies for equal masses and equal 
frequencies are consistently lower than those reported in 
Ref. [13. For A^ = 20, e.g., we find E = 41.35(8)to 
while Ref. [l^l reports 43.2(4)ftw. We speculate that this 
discrepancy can be traced back to the nodal structure 
of the trial wave function employed, and possibly also to 
the larger range of the two-body potential employed in 
Ref. 27 1 . In agreement with Ref. [131 , we find that the 
energies at unitarity show no shell structure. 

Using the LDA, which should be valid for sufficiently 
large A^, we relate the energy of the trapped Fermi system 
at unitarity to the universal parameter and E^j, i.e.. 



EiN) 



NI- 



(12) 



The parameter connects the energy per particle 
Ehom/N of the homogeneous system at unitarity and 
the energy per particle Epo of the homogeneous non- 
interacting Fermi gas, i.e., Ehom/N = ^kEfg- Here, we 
assumed that the functional dependence of Eiio„i/N on 
EpQ is the same for equal and unequal masses but that 
the universal parameter ^ depends on k. Applying the 
extended Thomas- Fermi (ETF) model to E^i [2g|, the 
energy of the trapped system at unitarity becomes 



EiN) 



(37V)4/3 



1 



(3A^) 



-2/3 



,(13) 



where = 1. The first term in Eq. (|13p is often referred 
to as Thomas-Fermi (TF) approximation. We also at- 
tempted to fit our FN-DMC energies by functional forms 
different from Eq. ()13|) . which included higher-order cor- 
rection terms or terms with other powers of A^; however, 
none of the alternative functional forms considered im- 
proved the description of our numerical results. Fitting 
our equal mass energies to Eq. treating as a pa- 
rameter, we find ^1 = 0.465. Our ^1 extracted from the 
trapped system is about 10% larger than that determined 
for the bulk system, i.e., ^1 = 0.42(1) [H,!!^, suggesting 
that one has to go to somewhat larger trapped systems 
to extrapolate the bulk ^1 with high accuracy within the 
LDA. Although the ^1 obtained from the fit to the ener- 
gies of the trapped system is larger than the correspond- 
ing bulk value, it is worthwhile noting that the simple 
functional form given in Eq. (|13p provides an excellent 
description of the energies of the trapped system. 

For K = 8, we find that our energies are best described 
if we treat and cg as fitting parameters, yielding = 
0.417 and cg = 0.27. The decrease of cg compared to the 
value predicted by the ETF model is most likely related 
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FIG. 6; (Color online) as a function of k. Inset: Ck as 
a function of k. The errorbars indicate the uncertainty of 
the fitting parameters; this uncertainty does not include the 
statistical errors of the FN-DMC energies. 



to the fact that the densities of the unequal mass species 
do, in contrast to the LDA treatment employed to derive 
Eq. (jl3p . not fully overlap at unitarity. Our calculations 
suggest that the bulk value for is somewhat smaller 
for K = 8 than for n = 1. To investigate this further, 
we additionally consider the energetics of systems with 
K = 4, 12, 16 and 20. Figure [5] shows the resulting 
(main panel) and (inset) extracted from our energies 
for A'^ = 2—20 fermions as a function of k. Both ^„ and 
vary smoothly, decrease with increasing k, and seem to 
approach a constant for large k. Furthermore, changes 
sign from positive to negative for k « 10. We emphasize 
that Eq. (fT3)) provides a rather good description of the 
energetics for all n thus empirically motivating the non- 
constant coefficient of the correction term. 

The fact that decreases with increasing k is in agree- 
ment with recent FN-DMC calculations for the homoge- 
neous system ^2^. However, this decrease is more pro- 
nounced for the trapped system than for the homoge- 
neous system. Standard BCS mean-field theory, in con- 
trast, predicts that the parameter ^, which determines 
many properties of dilute homogeneous two-com pon ent 
Fermi gases at unitarity , is independent of k |31| . 

We now comment further on the choice of the guid- 
ing function used to obtain the energies for 2 — 20 
that enter into our determination of and c^. As men- 
tioned earlier, trimer states with negative energy exist 
for K — 16 and 20. If we use the guiding function 'ipT2 
to model these systems, many-body configurations that 
contain three particles in close proximity are being sam- 
pled, giving rise to negative energies at unitarity. How- 
ever, if we use the guiding function ipri, i-e., if we con- 
struct the nodal surface by pairing spin-up and spin-down 
particles, many-body configurations that contain tightly- 
bound trimer states are not being sampled. Thus, the 
guiding function ipTi allows for a numerically stable char- 
acterization of a state with positive energy that has the 
same symmetry as ipri ■ We note that the decrease of 



with increasing k is already present in the energies for 
the two-body system for which we can determine the en- 
ergetics essentially exactly. This provides some evidence 
that the behaviors discussed in this section for unequal 
masses are not an artefact of our choice of guiding func- 
tion. 

The difference between the guiding functions ipri and 
'0T2 can also be understood from a different point of view. 
In the limit of vanishing confinement, the oscillator states 
used to construct ipT approach free-particle states. In 
this case, the nodal surface of the guiding function ipxi 
is compatible with a superfluid state, and the guiding 
function ^pT2 with a normal state [l^ [2^, [H, [s^l • Using 
this analogy, our results suggest that even fairly small 
trapped Fermi systems are better described by a "super- 
fluid wave function" than by a "normal wave function" . 
As A^ increases, the difference between the energies ob- 
tained for the superfluid and normal wave functions in- 
creases, presumably approaching the bulk values in the 
large A^ limit (Ehom/N = 0A2{1)Efg 0, HI for the 
superfluid state and Ehom/N = Q.^AEpc for the normal 
state 129,, 



IV. CONCLUSION 

This paper characterizes the BEC-BCS crossover 
physics of trapped two-component Fermi gases with 
varying mass ratio. Our results are obtained by solv- 
ing the stationary many-body Schrodinger equation for 
short-range model potential by two complementary ap- 
proaches. For the four-particle system, an essentially ex- 
act basis set expansion type technique, a CG approach, is 
used to determine the complete = 0^ spectrum. For 
up to A = 20 particles, the FN-DMC approach is used 
to determine upper bounds for energy of the BEC-BCS 
crossover branch. Treating the four-body system is chal- 
lenging, and interesting in its own right: The four-body 
system is the smallest non-trivial system exhibiting BEC- 
BCS crossover-like physics. Furthermore, the lessons 
learned from the four-body system aid the study of larger 
systems. Solving the Schrodinger equation for more than 
a few fermions by first principle methods is, despite the 
increasing available computer power, still a challenging 
task. In fact, it may be argued that Monte Carlo meth- 
ods are the only methods suitable. Unfortunately how- 
ever, assessing the accuracy of the assumptions going into 
Monte Carlo calculations, such as the nodal surface em- 
ployed in the FN-DMC approach, remains a challenge. 
Our calculations, which use the CG and FN-DMC ap- 
proaches in parallel, benchmark the strengths and lim- 
itations of the nodal surface employed in the FN-DMC 
calculations. 

From our four-body calculations in the deep BEG 
regime, we determine the scattering length add and effec- 
tive range rdd of the dimer-dimer system for two purely 
attractive short-range two-body potentials as a function 
of the mass ratio n. For up to k « 13.6, our dimer-dimer 
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scattering lengths add agree well with the values calcu- 
lated by Petrov et al. [1| . Our four-body calculations ex- 
tend beyond this mass ratio k and suggest that the ener- 
getically lowest-lying dimcr-dimer state varies smoothly 
as a function of k. We find that the energies of the dimer- 
dimer states for large k depend, just as for small k, at 
most weakly on the details of the two-body potential; 
we take this as numerical evidence that at least some 
properties of the dimer-dimer states are universal even if 
K, exceeds the value of 13.6. However, other properties 
of systems with large k such as the system's stability, 
encapsulated in the three-body recombination rate, are 
presumably controlled by the details of the short-range 
potential. These non-universal properties should be in- 
vestigated in the future. Also, future studies will have to 
investigate whether the comparatively large value of rdd 
can be measured indirectly by, e.g., a careful analysis of 
the density profile in the BEG regime. 

We also present calculations in the strongly-interacting 



unitarity regime for different mass ratios. Our calcula- 
tions for iV = 2 — 20 fermions show no shell structure. 
Application of the LDA to the trapped system implies 
that the universal parameter k depends weakly on the 
mass ratio. Our energies at unitarity for various mass 
ratios may aid in developing and refining numerically less 
demanding treatments of two-component Fermi gases. 
Recently, e.g., Bulgac [9,] proposed a density functional 
theory applicable to trapped equal-mass two-component 
Fermi gases at unitarity. Our results presented here pro- 
vide much needed benchmarks for such theories. Our un- 
equal mass studies present the first first principle treat- 
ment of such systems under confinement. Our analysis 
provides a first step towards a deeper understanding of 
these systems, but much room for further investigations, 
including the investigation of connections to mean-field 
treatments [H [M HI, [H, [13, [H, ill , remains. 

We gratefully acknowledge discussions with S. Giorgini 
and S. Rittenhouse, and support by the NSF. 
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